###############################################################################
###   Raz (2024).                 													                ###
###   Soil Heterogeneity, Social Learning, and the Formation of Close-knit  ###
###   Communities.                                                          ###
###-------------------------------------------------------------------------###
###   Make Tables       								                                    ###
###############################################################################

# Author: Itzchak Tzachi Raz
# Last Revised: 11/30/2024

# Note: some input data files could not be made publicly available.  


# Preambles ---------------------------------------------------------------

#---------------------------------
# --- Load packages  
#---------------------------------

library(dplyr)
library(stringr)
library(data.table)
library(lfe)
library(stargazer)
library(robomit)
library(compareGroups)
library(expss)

#---------------------------------
# --- Paths
#---------------------------------

scriptPath <- dirname(rstudioapi::getSourceEditorContext()$path) 
if(dir.exists(scriptPath) == FALSE) {
  print("ERROR in script path")
}
BasicPath <- dirname(scriptPath)
DataPath <- paste0(BasicPath, "/data")
PaperTablesPath <- paste0(BasicPath, "/output/paper/tables")
AppendixTablesPath <- paste0(BasicPath, "/output/appendix/tables")

#---------------------------------
# -- Parameters
#---------------------------------

my.GeoCtrls <- "Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Shoreline + dist2Lakes"
my.SmoothLocCtrls <- "Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Shoreline + dist2Lakes + XCoord + YCoord + XCoord2 + YCoord2 + XYCoord"
my.AgSuitabilityCtrls <- "Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Shoreline + dist2Lakes + XCoord + YCoord + XCoord2 + YCoord2 + XYCoord + AlfalfaSI + CottonSI + MaizeSI + OatSI + RyeSI + SugercaneSI + Sweet_potatoSI + TobaccoSI +  WheatSI + White_potatoSI"
my.HigherOrderCtrls <- "Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Shoreline + dist2Lakes + XCoord + YCoord + XCoord2 + YCoord2 + XYCoord + AlfalfaSI + CottonSI +  MaizeSI + OatSI + RyeSI + SugercaneSI + Sweet_potatoSI + TobaccoSI +  WheatSI + White_potatoSI + ElevationSD + SlopeSD + Flow_accSD + PrecipitationSD + TemperatureSD + AlfalfaSI_SD + CottonSI_SD + MaizeSI_SD + OatSI_SD + RyeSI_SD + SugercaneSI_SD + Sweet_potatoSI_SD + TobaccoSI_SD +  WheatSI_SD + White_potatoSI_SD"

my.colClasses <- c("STATE_TERR"="factor", "YEAR"="factor")


# -------------------------------------------------------------------------
# ----------                    MAIN TABLES                      ----------
# -------------------------------------------------------------------------

# Table 1 -----------------------------------------------------------------

table1 <- function(dep.var.name, indep.var.name, FEs, reg.data, clusterVar,
                   keep.vars = "SHI25", OsterDelta = FALSE) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " | 1 | 0 | ", clusterVar)), 
             data = reg.data)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " | ", FEs, " | 0 | ", clusterVar)), 
             data = reg.data)
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.GeoCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.AgSuitabilityCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.HigherOrderCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  if(OsterDelta) {
    reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                           float = FALSE,
                           header = FALSE,
                           model.names = FALSE,
                           model.numbers = FALSE,
                           keep = keep.vars,
                           keep.stat = c("n", "rsq"),
                           add.lines = list(c("Oster $\\delta$ for $\\beta=0$", "", 
                                              o_delta(y = dep.var.name,
                                                      x = "SHI25",
                                                      con = paste0("+ factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m6)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2), 
                                              o_delta(y = dep.var.name,
                                                      x = "SHI25",
                                                      con = paste0("+ ", my.GeoCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m6)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2), 
                                              o_delta(y = dep.var.name,
                                                      x = "SHI25",
                                                      con = paste0("+ ", my.SmoothLocCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m6)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2), 
                                              o_delta(y = dep.var.name,
                                                      x = "SHI25",
                                                      con = paste0("+ ", my.AgSuitabilityCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m6)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2), 
                                              o_delta(y = dep.var.name,
                                                      x = "SHI25",
                                                      con = paste0("+ ", my.HigherOrderCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m6)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2))),
                           table.layout = "tas")
    
    return(reg.table)
  } else{
    reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                           float = FALSE,
                           header = FALSE,
                           model.names = FALSE,
                           model.numbers = FALSE,
                           keep = keep.vars,
                           keep.stat = c("n", "rsq"),
                           table.layout = "tas")
    
    return(reg.table)
  }
}

FilePath <- paste0(DataPath, "/final/CountyLevelData.csv")
CountyLevelData <- fread(FilePath, colClasses = my.colClasses)


#--- Panel A: LNI ---
reg_out <- table1(dep.var.name = "LNI_county_Native",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:8])

outfile <- paste0(PaperTablesPath, "/table_1A.tex")
write(reg_out, outfile)

#--- Panel B: ICM ---
reg_out <- table1(dep.var.name = "ICM_Native",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:8])

outfile <- paste0(PaperTablesPath, "/table_1B.tex")
write(reg_out, outfile)

#--- Panel C: TNI ---
reg_out <- table1(dep.var.name = "TNI_Native",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:8])

outfile <- paste0(PaperTablesPath, "/table_1C.tex")
write(reg_out, outfile)


#--- Panel D: RHI ---
reg_out <- table1(dep.var.name = "RHI",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:8])

outfile <- paste0(PaperTablesPath, "/table_1D.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
  c(
  paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
  paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
  paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
  paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
  paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
  collapse = "\n")  

outfile <- paste0(PaperTablesPath, "/table_1_foot.tex")
write(footLines, outfile)

# Table 2 -----------------------------------------------------------------
# Selective-in migration

table2 <- function(dep.var.name, indep.var.name, FEs, ind.ctrls, reg.data) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " | 1 | 0 | GISJOIN")), 
             data = reg.data)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, 
                            " | ", FEs, " | 0 | GISJOIN")), 
             data = reg.data)
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.GeoCtrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.AgSuitabilityCtrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.HigherOrderCtrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m7 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.HigherOrderCtrls,
                            " | ", FEs, " + ", ind.ctrls, " | 0 | GISJOIN")),
             data = reg.data)
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  return(c(reg.table[3:4], "\\\\[-1.8ex] ", reg.table[6:7]))
}

#--- Panel A: farming experience ---
FilePath <- paste0(DataPath, "/final/SelectionFarmingData.csv")
SelectionFarmingData <-fread(FilePath, colClasses = c("STATE_TERR_y"="factor",
                                                      "STATE_TERR_z"="factor",
                                                      "YEAR"="factor",
                                                      "GISJOIN_z" = "factor"))

reg_out <- table2(dep.var.name = "farmer_z",
                  indep.var.name = "SHI25",
                  FEs = "GISJOIN_z:STATE_TERR_y:YEAR",
                  ind.ctrls = "age_y",
                  reg.data = SelectionFarmingData)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

outfile <- paste0(PaperTablesPath, "/table_2A.tex")
write(reg_out, outfile)

#--- Panel B: LNI ---

FilePath <- paste0(DataPath, "/final/migrants/MigrantsNative_2.csv")
MigrantsData <- fread(FilePath, 
                      colClasses = c("statefips"="factor", "YEAR"="factor",
                                     "sex" ="factor", "birthYear"="factor",
                                     "birthOrder"="factor",
                                     "b_statefips" = "factor"))

reg_out <- table2(dep.var.name = "LNI_state_Native",
                  indep.var.name = "SHI25",
                  FEs = "b_statefips:statefips:YEAR",
                  ind.ctrls = "sex:birthYear + birthOrder",
                  reg.data = MigrantsData[Post==0])

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

outfile <- paste0(PaperTablesPath, "/table_2B.tex")
write(reg_out, outfile)


#--- Panel C: ICM ---

reg_out <- table2(dep.var.name = "ICM",
                  indep.var.name = "SHI25",
                  FEs = "b_statefips:statefips:YEAR",
                  ind.ctrls = "sex:birthYear + f_age",
                  reg.data = MigrantsData[, .SD[1], by = 'familyFE'])

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

outfile <- paste0(PaperTablesPath, "/table_2C.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Origin $\\times$ Destination State", paste(rep(" ", 8), collapse = " & "), "\\\\"),
      paste(" \\; \\; \\; \\; \\; \\  $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"), 
      paste("Individual Controls", paste(rep(" ", 8), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(PaperTablesPath, "/table_2_foot.tex")
write(unlist(footLines), outfile)

# Table 3 -----------------------------------------------------------------
# Migrants, DD and DDD

table3 <- function(reg.data, ind.ctrls = NA) {
  
  if(is.na(ind.ctrls)) {
    m1 <- felm(LNI_state_Native ~ Post:SHI25 | familyFE + Rtime | 0 | GISJOIN, 
               data = reg.data)
    
    m2 <- felm(LNI_state_Native ~ Post:SHI25 | familyFE + Rtime | 0 | GISJOIN, 
               data = reg.data[farmerDad==1])
    
    m3 <- felm(LNI_state_Native ~ Post:SHI25 | familyFE + Rtime | 0 | GISJOIN, 
               data = reg.data[farmerDad==0])
    
    m4 <- felm(LNI_state_Native ~ Post:farmerDad:SHI25 + Post:farmerDad + Post:SHI25 | 
                 familyFE + Rtime | 0 | GISJOIN, 
               data = reg.data)
    
  } else {
    m1 <- felm(formula(paste0("LNI_state_Native ~ Post:SHI25 | familyFE + Rtime + ",
                              ind.ctrls, " | 0 | GISJOIN")), 
               data = reg.data)
    
    m2 <- felm(formula(paste0("LNI_state_Native ~ Post:SHI25 | familyFE + Rtime  + ",
                              ind.ctrls, "| 0 | GISJOIN")), 
               data = reg.data[farmerDad==1])
    
    m3 <- felm(formula(paste0("LNI_state_Native ~ Post:SHI25 | familyFE + Rtime  + ",
                              ind.ctrls, "| 0 | GISJOIN")), 
               data = reg.data[farmerDad==0])
    
    m4 <- felm(formula(paste0("LNI_state_Native ~ Post:farmerDad:SHI25 + Post:farmerDad + Post:SHI25 | 
                 familyFE + Rtime  + ",ind.ctrls, "| 0 | GISJOIN")), 
               data = reg.data)
  }

  reg.table <- stargazer(m1,m2,m3,m4,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = c("Post:SHI25", "Post:farmerDad:SHI25"),
                         keep.stat = c("n", "rsq"),
                         add.lines = c(list(c("Dependent Variable Mean", 
                                              mean(m1$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m2$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m3$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m4$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ",")),
                                            c("Dependent Variable SD", 
                                              sd(m1$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m2$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m3$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m4$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ",")))),
                         table.layout = "tas")

  reg.table <- str_replace(reg.table, "Post:SHI25", 
                           "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\Post Migration  $\\\\times$ \\\\\\\\ \\\\quad Soil Heterogeneity}}")
  reg.table <- str_replace(reg.table, "Post:farmerDad:SHI25", 
                           "\\\\multirow[t]{3}{*}{\\\\makecell[l]{\\\\\\\\ \\\\\\\\Post Migration  $\\\\times$ \\\\\\\\ \\\\quad  Farmers$'$ Household  $\\\\times$ \\\\\\\\ \\\\qquad Soil Heterogeneity}}")

  
  return(c(reg.table[3:8], "& & & & \\\\ \\hline \\\\[-1.8ex]  \\\\[-1.8ex] ", reg.table[9:12]))
}


reg_out <- table3(reg.data = MigrantsData)
outfile <- paste0(PaperTablesPath, "/table_3.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Family Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Relative YOB Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(PaperTablesPath, "/table_3_foot.tex")
write(unlist(footLines), outfile)

# Table 4 -----------------------------------------------------------------
# Long-run impact on migrants' children

table4 <- function(dep.var.name, FEs, reg.ctrls, reg.data) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", reg.ctrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", reg.ctrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data[farmerDad==1])
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", reg.ctrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data[farmerDad==0])
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ SHI25*farmerDad + ", reg.ctrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  reg.table <- stargazer(m1,m2,m3,m4,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         add.lines = c(list(c("Dependent Variable Mean", 
                                              mean(m1$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m2$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m3$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m4$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ",")),
                                            c("Dependent Variable SD", 
                                              sd(m1$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m2$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m3$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m4$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ",")))),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "SHI25:farmerDad", 
                           "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\Soil Heterogeneity  $\\\\times$ \\\\\\\\ \\\\quad Farmers$'$ Household}}")
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  
  return(c(reg.table[3:8], "\\\\[-1.8ex] ", reg.table[9:12]))
}

FilePath <- paste0(DataPath, "/final/MigrantsChildrenData.csv")
MigrantsChildren <- fread(FilePath)

MigrantsChildren[, ':=' (YEAR = as.factor(YEAR), 
                         age = as.factor(age), 
                         statefip_y = as.factor(statefip_y),
                         bpl = as.factor(bpl),
                         ICM = as.factor(ICM),
                         highLNI = as.factor(highLNI), 
                         Rtime = as.factor(Rtime))]


#--- Panel A: Remained in the destination county ---
reg_out <- table4(dep.var.name = "percent",
                  FEs = "bpl:statefip_y:YEAR:Rtime:highLNI:ICM",
                  reg.ctrls = my.SmoothLocCtrls,
                  reg.data = MigrantsChildren)

outfile <- paste0(PaperTablesPath, "/table_4A.tex")
write(reg_out, outfile)


# --- Panel B: Married a spouse born in the destination state ---
reg_out <- table4(dep.var.name = "ICM_d",
                  FEs = "bpl:statefip_y:YEAR:Rtime:highLNI:ICM",
                  reg.ctrls = my.SmoothLocCtrls,
                  reg.data = MigrantsChildren)

outfile <- paste0(PaperTablesPath, "/table_4B.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Cell Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(PaperTablesPath, "/table_4_foot.tex")
write(unlist(footLines), outfile)

# Table 5 -----------------------------------------------------------------
# SHI and farmers learning

table5 <- function(dep.var.name, learning.var, high.learning.var, FEs, 
                   reg.ctrls, reg.data) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", reg.ctrls,
                            " | ", FEs, " | 0 | grid100m")),
             data = reg.data)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ ", high.learning.var, " + ", reg.ctrls,
                            " | ", FEs, " | 0 | grid100m")),
             data = reg.data)
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", high.learning.var, " + ", reg.ctrls,
                            " | ", FEs, " | 0 | grid100m")),
             data = reg.data)
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", reg.ctrls,
                            " | ", FEs, " | 0 | grid100m")),
             data = reg.data %>% filter(eval(parse(text = paste0(high.learning.var, "== 1")))))
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", reg.ctrls,
                            " | ", FEs, " | 0 | grid100m")),
             data = reg.data %>% filter(eval(parse(text = paste0(high.learning.var, "== 0")))))
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ SHI25*", high.learning.var, " + ", reg.ctrls,
                            " | ", FEs, " | 0 | grid100m")),
             data = reg.data)
  
  m7 <- felm(formula(paste0(learning.var, " ~ SHI25 + ", reg.ctrls,
                            " | ", FEs, " | 0 | grid100m")),
             data = reg.data[YEAR==1930])
  
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = c("SHI25", high.learning.var, paste0("SHI25:",high.learning.var)),
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, paste0("SHI25:",high.learning.var), 
                           "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\Soil Heterogeneity  $\\\\times$ \\\\\\\\ \\\\quad High Learning Potential}}")
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  reg.table <- str_replace(reg.table, high.learning.var, "High Learning Potential")
  
  return(c(reg.table[3:11], "\\\\[-1.8ex] ", reg.table[12:13]))
}


#--- Panel A: Growth in Fertilizers Use ---
FilePath <- paste0(DataPath, "/final/FertilizerData.csv")
FertilizerData <- fread(FilePath, colClasses = my.colClasses)

reg_out <- table5(dep.var.name = "g_shareFertilizer",
                  FEs = "STATE_TERR:YEAR",
                  learning.var = "gains5", 
                  high.learning.var = "high.gains5",
                  reg.ctrls = my.SmoothLocCtrls,
                  reg.data = FertilizerData)

outfile <- paste0(PaperTablesPath, "/table_5A.tex")
write(reg_out, outfile)

#--- Panel B: Growth in Wheat Share ---
FilePath <- paste0(DataPath, "/final/WheatShareData.csv")
WheatShareData <- fread(FilePath, colClasses = my.colClasses)

reg_out <- table5(dep.var.name = "g_shareWheat",
                  FEs = "STATE_TERR:YEAR",
                  learning.var = "WheatDiff", 
                  high.learning.var = "high.gains",
                  reg.ctrls = my.SmoothLocCtrls,
                  reg.data = WheatShareData)

outfile <- paste0(PaperTablesPath, "/table_5B.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(PaperTablesPath, "/table_5_foot.tex")
write(unlist(footLines), outfile)

# Table 6 ---------------------------------------------------------------
# Selective Out-Migration

FilePath <- paste0(DataPath, "/final/OutMigrationData.csv")
OutMigrationData <- fread(FilePath, colClasses = my.colClasses)

m1 <- felm(formula(paste0("percent ~ SHI25 + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData)

m2 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData)

m3 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData[highLNI==0])

m4 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData[highLNI==1])

m5 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData[ICM==0])

m6 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData[ICM==1])

m7 <- felm(formula(paste0("percent ~ SHI25*farmer_y*highLNI + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData)


reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                       float = FALSE,
                       header = FALSE,
                       model.names = FALSE,
                       model.numbers = FALSE,
                       keep = c("^SHI25$", "^farmer_y$", "^SHI25:farmer_y$", "^SHI25:farmer_y:highLNI$"),
                       keep.stat = c("n", "rsq"),
                       table.layout = "tas")

reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
reg.table <- str_replace(reg.table, "farmer\\\\_y", "Farmer")
reg.table <- str_replace(reg.table, "Soil Heterogeneity:Farmer:highLNI", 
                         "\\\\multirow[t]{4}{*}{\\\\makecell[l]{\\\\\\\\ \\\\\\\\ \\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer  $\\\\times$ \\\\\\\\ \\\\qquad High Communal \\\\\\\\ \\\\qquad Identification}}")
reg.table <- str_replace(reg.table, "Soil Heterogeneity:Farmer", 
                         "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer}}")


outfile <- paste0(PaperTablesPath, "/table_6.tex")
write(c(reg.table[3:14],
        "  & & & & & & \\\\ \\\\", "\\hline \\\\[-1.8ex]  \\\\[-1.8ex] ",
        reg.table[15:16]), outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(PaperTablesPath, "/table_6_foot.tex")
write(unlist(footLines), outfile)


# -------------------------------------------------------------------------
# ----------                  APPENDIX TABLES                    ----------
# -------------------------------------------------------------------------


# Appendix Table A.1 ------------------------------------------------------
# SHI and geo-climatic covariates

# Adds starts for significance level in regression tables
my.stars <- function(reg = NA, reg.num = NA) {
  
  stars.text <- ""
  if(abs(reg$beta[reg.num]/reg$cse[reg.num]) > 2.326) {
    stars.text <- "$^{***}$"
  } else if(abs(reg$beta[reg.num]/reg$cse[reg.num]) > 1.96) {
    stars.text <- "$^{**}$"
  } else if(abs(reg$beta[reg.num]/reg$cse[reg.num]) > 1.645) {
    stars.text <- "$^{*}$"
  }
  return(stars.text)
}

# Run regression w and w/o controls and saves output for table
tableA1.line.short <- function(my.covar = NA, num.base = NA, covar.label = NA) {
  
  reg_base <- felm(paste0("SHI25 ~  ", my.covar, " | 1 | 0 | grid100m") %>%
                     formula(), data = geoData) 
  reg_state <- felm(paste0("SHI25 ~  ", my.covar, " | STATE_TERR | 0 | grid100m") %>%
                      formula(), data = geoData) 
  reg_smooth <- felm(paste0("SHI25 ~  ", my.covar, " +
                             XCoord + YCoord + XCoord2 + YCoord2 + XYCoord | 
                             1 | 0 | grid100m") %>%
                       formula(), data = geoData) 
  reg_smooth_state <- felm(paste0("SHI25 ~  ", my.covar, " +
                                   XCoord + YCoord + XCoord2 + YCoord2 + XYCoord | 
                                   STATE_TERR | 0 | grid100m") %>%
                             formula(), data = geoData) 
  
  
  reg.out <- c(paste0("\\\\[-1.8ex] ", covar.label, " & ", 
                      reg_base$beta[2] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","), 
                      my.stars(reg_base, 2),
                      " & ", 
                      reg_state$beta[1] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","),
                      my.stars(reg_state, 1),
                      " & ", 
                      reg_smooth$beta[2] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","),
                      my.stars(reg_smooth, 2),
                      " & ", 
                      
                      all$beta[num.base+1] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","), 
                      my.stars(all, num.base+1),
                      " & ", 
                      all_smooth_state$beta[num.base] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","),
                      my.stars(all_smooth_state, num.base),
                      " \\\\"),
               
               paste0(" & (", 
                      reg_base$cse[2] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","), 
                      ") & (",  
                      reg_state$cse[1] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","), 
                      ") & (",
                      reg_smooth$cse[2] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","), 
                      ") & (",  
                      
                      all$cse[num.base+1] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","), 
                      ") & (",  
                      
                      all_smooth_state$cse[num.base] %>% 
                        formatC(format="fg", digits = 2, big.mark = ","),
                      ") \\\\"))
  
  return(reg.out)
}


FilePath <- paste0(DataPath, "/final/LongRunData.csv")
geoData <- fread(FilePath)
keepVarList <- c("GISJOIN", "STATE_TERR", "SHI25", "grid100m",
                 "Elevation", "Slope", "Flow_acc", "Precipitation", "Temperature", 
                 "RiverDensity", "maxval", "SHAPE_AREA", "dist2NavRiver", "dist2Lakes", 
                 "dist2Shoreline", "XCoord", "YCoord", "XCoord2", "YCoord2", "XYCoord", 
                 "AlfalfaSI", "CottonSI",  "MaizeSI", "OatSI", "RyeSI", "SugercaneSI", 
                 "Sweet_potatoSI", "TobaccoSI",  "WheatSI", "White_potatoSI",
                 "ElevationSD", "SlopeSD", "Flow_accSD", "PrecipitationSD", "TemperatureSD",
                 "AlfalfaSI_SD", "CottonSI_SD", "MaizeSI_SD", "OatSI_SD", "RyeSI_SD", 
                 "SugercaneSI_SD", "Sweet_potatoSI_SD", "TobaccoSI_SD", "WheatSI_SD", 
                 "White_potatoSI_SD")
geoData <- geoData[, .SD, .SDcols = keepVarList]


geoData <- geoData %>%
  mutate(Elevation = scale(Elevation),
         Slope = scale(Slope),
         Flow_acc = scale(Flow_acc),
         Precipitation = scale(Precipitation),
         Temperature = scale(Temperature),
         RiverDensity = scale(RiverDensity),
         maxval = scale(maxval),
         SHAPE_AREA = scale(SHAPE_AREA),
         dist2NavRiver = scale(dist2NavRiver),
         dist2Lakes = scale(dist2Lakes),
         dist2Shoreline = scale(dist2Shoreline))


# Run regression with full geo controls
all <- felm(SHI25 ~ Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Lakes + dist2Shoreline + dist2Lakes| 
              1 | 0 | grid100m, data = geoData) 
all_state <- felm(SHI25 ~ Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Lakes  + dist2Shoreline + dist2Lakes | 
                    STATE_TERR | 0 | grid100m, data = geoData) 
all_smooth <- felm(SHI25 ~ Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Lakes  + dist2Shoreline + dist2Lakes +
                     XCoord + YCoord + XCoord2 + YCoord2 + XYCoord | 
                     1 | 0 | grid100m, data = geoData) 
all_smooth_state <- felm(SHI25 ~  
                           Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Lakes  + dist2Shoreline + dist2Lakes +
                           XCoord + YCoord + XCoord2 + YCoord2 + XYCoord | 
                           STATE_TERR | 0 | grid100m, data = geoData) 

reg.out <- c(
  tableA1.line.short("Temperature", 5, "Temperature"),
  tableA1.line.short("Precipitation", 4, "Precipitation"),
  tableA1.line.short("Slope", 2, "Slope"),
  tableA1.line.short("Elevation", 1, "Elevation"),
  tableA1.line.short("maxval", 7, "Productivity"),
  tableA1.line.short("Flow_acc", 3, "Flow accumulation"),
  tableA1.line.short("RiverDensity", 6, "River density"),
  tableA1.line.short("SHAPE_AREA", 8, "Total area"),
  tableA1.line.short("dist2NavRiver", 9, "Distance to navigated rivers"),
  tableA1.line.short("dist2Lakes", 9, "Distance to lakes"),
  tableA1.line.short("dist2Shoreline", 10, "Distance to shoreline"),
  " & & & & &  \\\\",
  "\\hline \\\\[-1.8ex] \\\\[-1.8ex] ",
  "State Fixed Effects      &  & \\checkmark &  & & \\checkmark \\\\",
  "Smooth Location Controls &  &  & \\checkmark & & \\checkmark \\\\",
  "Geoclimatic Controls     &  &  &  & \\checkmark & \\checkmark \\\\" 
)

outfile <- paste0(AppendixTablesPath, "/table_A_1.tex")
write(reg.out, outfile)

rm(all, all_state, all_smooth, all_smooth_state)


# Appendix Table A.2 ------------------------------------------------------
# Validation of the SHI

#--- Panel A: Variation of Agriculture Suitability ---
AgSuitabilityISs <- c("AlfalfaSI_SD", "CottonSI_SD", "MaizeSI_SD", "OatSI_SD", "RyeSI_SD", "SugercaneSI_SD", "Sweet_potatoSI_SD", "TobaccoSI_SD", "WheatSI_SD", "White_potatoSI_SD")

for (SI in AgSuitabilityISs) {
  geoData[[(paste0("z", SI))]] <- scale(eval(parse(text = paste0("geoData$", SI)))) %>% as.numeric()
}

AgSuitabilityISs <- paste0("z", AgSuitabilityISs)
geoData[, GAEZheterogeneity := sum(.SD, na.rm=T)/length(AgSuitabilityISs), 
         by = "GISJOIN", .SDcols = AgSuitabilityISs]

my.HigherOrderCtrlsTEMP <- my.HigherOrderCtrls
my.HigherOrderCtrls <- "Elevation + Slope + Flow_acc + Precipitation + Temperature + RiverDensity + maxval + SHAPE_AREA + dist2NavRiver + dist2Shoreline + dist2Lakes + XCoord + YCoord + XCoord2 + YCoord2 + XYCoord + AlfalfaSI + CottonSI +  MaizeSI + OatSI + RyeSI + SugercaneSI + Sweet_potatoSI + TobaccoSI +  WheatSI + White_potatoSI + ElevationSD + SlopeSD + Flow_accSD + PrecipitationSD + TemperatureSD"

reg_out <- table1(dep.var.name = "GAEZheterogeneity",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR",
                  reg.data = geoData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:5], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_A_2A.tex")
write(reg_out, outfile)

my.HigherOrderCtrls <- my.HigherOrderCtrlsTEMP
rm(my.HigherOrderCtrlsTEMP, geoData)


#--- Panel B: Agricultural diversity ---
FilePath <- paste0(DataPath, "/final/AgriculturalDiversityData.csv")
AgDivData <- fread(FilePath, colClasses = my.colClasses)

reg_out <- table1(dep.var.name = "AgDiversity",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = AgDivData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:5], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_A_2B.tex")
write(reg_out, outfile)

rm(AgDivData)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_A_2_foot.tex")
write(unlist(footLines), outfile)



# Appendix Table B.1 ------------------------------------------------------

FilePath <- paste0(DataPath, "/final/LNIExampleData.csv")
LNIexampleData <- fread(FilePath)

my.table <- c(
  paste0("\\begin{tabular}{@{\\extracolsep{20pt}}lcccc}"),
  paste0("\\\\[-1.8ex]\\hline"),
  paste0("\\hline \\\\[-1.8ex]"), 
  paste0("\\\\[-1.8ex] & \\multicolumn{2}{c}{Arkansas} & \\multicolumn{2}{c}{Massachusetts} \\\\"),
  paste0("\\cline{2-3} \\cline{4-5}"),
  paste0("\\\\[-1.8ex]"), 
  paste0("\\\\[-1.8ex]  & Share &  LNI & Share &  LNI \\\\"),
  paste0("\\cline{2-2} \\cline{3-3} \\cline{4-4} \\cline{5-5}"),
  paste0("\\\\[-1.8ex]"),
  
  paste0("\\\\[-1.8ex] ", " & ", "\\multicolumn{4}{c}{\\textit{Panel A: Common in one, rare in the other}} \\\\"), 
  paste0("\\cline{2-5} \\\\[-1.8ex]"), 
  
  paste0("\\\\[-1.8ex] Billie", " & ",
         LNIexampleData[first_name=="billie" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="billie" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="billie" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="billie" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] Floyd", " & ",
         LNIexampleData[first_name=="floyd" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="floyd" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="floyd" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="floyd" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] Jerry", " & ",
         LNIexampleData[first_name=="jerry" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="jerry" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="jerry" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="jerry" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  
  paste0("\\\\[-1.8ex] Francis",  " & ",
         LNIexampleData[first_name=="francis" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="francis" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="francis" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="francis" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] Peter",  " & ",
         LNIexampleData[first_name=="peter" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="peter" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="peter" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="peter" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] Frederick",  " & ",
         LNIexampleData[first_name=="frederick" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="frederick" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="frederick" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="frederick" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] \\\\[-1.8ex] ", " & ", "\\multicolumn{4}{c}{\\textit{Panel B: Common in both}} \\\\"), 
  paste0("\\cline{2-5} \\\\[-1.8ex]"), 
  
  paste0("\\\\[-1.8ex] James",  " & ",
         LNIexampleData[first_name=="james" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="james" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="james" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="james" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] John",  " & ",
         LNIexampleData[first_name=="john" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="john" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="john" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="john" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] William",  " & ",
         LNIexampleData[first_name=="william" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="william" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="william" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="william" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] Robert",  " & ",
         LNIexampleData[first_name=="robert" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="robert" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="robert" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="robert" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] Donald",  " & ",
         LNIexampleData[first_name=="donald" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="donald" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="donald" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="donald" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\\\[-1.8ex] \\\\[-1.8ex] ", " & ", "\\multicolumn{4}{c}{\\textit{Panel C: Rare in both}} \\\\"), 
  paste0("\\cline{2-5} \\\\[-1.8ex]"),   
  
  paste0("\\\\[-1.8ex] Walker",  " & ",
         LNIexampleData[first_name=="walker" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="walker" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="walker" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="walker" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2),
         " \\\\"),
  
  paste0("\\\\[-1.8ex] Wyatt",  " & ",
         LNIexampleData[first_name=="wyatt" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="wyatt" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         "--- & --- ",
         " \\\\"),
  
  paste0("\\\\[-1.8ex] Carleton",  " & ",
         LNIexampleData[first_name=="carleton" & statefips=="5", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="carleton" & statefips=="5", LNI] %>% 
           formatC(format="f", digits = 2), " & ",
         LNIexampleData[first_name=="carleton" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="carleton" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  
  paste0("\\\\[-1.8ex] Kevin",  " & ",
         "--- & --- & ",
         LNIexampleData[first_name=="kevin" & statefips=="25", nameShare] %>% 
           formatC(format="f", digits = 4), " & ",
         LNIexampleData[first_name=="kevin" & statefips=="25", LNI] %>% 
           formatC(format="f", digits = 2), 
         " \\\\"),
  
  paste0("\\hline"),
  paste0("\\hline \\\\[-1.8ex] "),
  
  paste0("\\end{tabular} "))

outfile <- paste0(AppendixTablesPath, "/table_B_1.tex")
write(my.table, outfile)

rm(LNIexampleData, my.table)

# Appendix Table C.1 ------------------------------------------------------
# Kinship Tightness

#--- Panel A: Median Kin Propinquity ---
reg_out <- table1(dep.var.name = "MKP_Native",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:8])

outfile <- paste0(AppendixTablesPath, "/table_C_1A.tex")
write(reg_out, outfile)

#--- Panel B: Kin Propinquity Rate ---
reg_out <- table1(dep.var.name = "KPR_Native",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:8])

outfile <- paste0(AppendixTablesPath, "/table_C_1B.tex")
write(reg_out, outfile)

#--- Panel C: Strength of Family Ties Index ---
reg_out <- table1(dep.var.name = "SFTI_Native",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:8])

outfile <- paste0(AppendixTablesPath, "/table_C_1C.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_C_1_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table C.2 ------------------------------------------------------
# LNI, Year-by-year

yearsList <- setdiff(seq.int(from = 1850, to = 1940, by = 10), 1890)

i <- 1
for (y in yearsList) {
  assign(paste0("m",i), felm(formula(paste0("LNI_county_Native ~ SHI25 + ", my.SmoothLocCtrls,
                                          " | STATE_TERR | 0 | grid100m")),
                           data = CountyLevelData[YEAR==y]))
  i <- i + 1
}

reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,m8,m9,
                       float = FALSE,
                       header = FALSE,
                       model.names = FALSE,
                       model.numbers = FALSE,
                       keep = "SHI25",
                       keep.stat = c("n", "rsq"),
                       add.lines = list(c("Dependent Variable Mean",  
                                          mean(m1$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ","), 
                                          mean(m2$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ","), 
                                          mean(m3$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ","), 
                                          mean(m4$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ","), 
                                          mean(m5$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ","),
                                          mean(m6$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ","),
                                          mean(m7$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ","),
                                          mean(m8$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ","),
                                           mean(m9$response) %>% 
                                                   formatC(format="f", digits = 2, big.mark = ",")),
                                        c("Dependent Variable SD",  
                                          sd(m1$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","), 
                                          sd(m2$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","), 
                                          sd(m3$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","), 
                                          sd(m4$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","), 
                                          sd(m5$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","),
                                          sd(m6$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","),
                                          sd(m7$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","),
                                          sd(m8$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","),
                                          sd(m9$response) %>% 
                                            formatC(format="f", digits = 2, big.mark = ","))),
                       table.layout = "tas")

reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")


outfile <- paste0(AppendixTablesPath, "/table_C_2.tex")
write(c(reg.table[3:5], "\\hline \\\\[-1.8ex]  \\\\[-1.8ex] ", reg.table[6:9]), 
      outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 9), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 9), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 9), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_C_2_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table C.3 ---------------------------------------------------------------
# Long-Run Features of Close-Knit Communities

FilePath <- paste0(DataPath, "/final/LongRunData.csv")
LongRunData <- fread(FilePath)

#--- Panel A: Social Clustering ---

reg_out <- table1(dep.var.name = "z_clustering_county",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR",
                  reg.data = LongRunData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_C_3A.tex")
write(reg_out, outfile)



#--- Panel B: Social Support Ratio ---
reg_out <- table1(dep.var.name = "z_support_ratio_county",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR",
                  reg.data = LongRunData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_C_3B.tex")
write(reg_out, outfile)


#--- Panel C: Communal Moral Values ---
FilePath <- paste0(DataPath, "/final/MFQData.csv")
MFQData <- fread(FilePath, colClasses = my.colClasses)

m1 <- felm(CMV ~ SHI25 | 
             gender + race + byear | 0 | grid100m, 
           weights = MFQData$res_ratio, 
           data = MFQData)

m2 <- felm(CMV ~ SHI25 | 
             gender + race + byear + STATE_TERR | 0 | grid100m, 
           weights = MFQData$res_ratio, 
           data = MFQData)

m3 <- felm(formula(paste0("CMV ~ SHI25 + ", my.GeoCtrls,
                          " | gender + race + byear + STATE_TERR | 0 | grid100m")),
           weights = MFQData$res_ratio, 
           data = MFQData)

m4 <- felm(formula(paste0("CMV ~ SHI25 + ", my.SmoothLocCtrls,
                          " | gender + race + byear + STATE_TERR | 0 | grid100m")),
           weights = MFQData$res_ratio, 
           data = MFQData)

m5 <- felm(formula(paste0("CMV ~ SHI25 + ", my.AgSuitabilityCtrls,
                          " | gender + race + byear + STATE_TERR | 0 | grid100m")),
           weights = MFQData$res_ratio, 
           data = MFQData)

m6 <- felm(formula(paste0("CMV ~ SHI25 + ", my.HigherOrderCtrls,
                          " | gender + race + byear + STATE_TERR | 0 | grid100m")),
           weights = MFQData$res_ratio, 
           data = MFQData)

reg_out <- stargazer(m1,m2,m3,m4,m5,m6,
                     float = FALSE,
                     header = FALSE,
                     model.names = FALSE,
                     model.numbers = FALSE,
                     keep = "SHI25",
                     keep.stat = c("n", "rsq"),
                     table.layout = "tas")

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])  

outfile <- paste0(AppendixTablesPath, "/table_C_3C.tex")
write(reg_out, outfile)


#--- Panel D: Excess Support for Trump in 2016 (Enke, 2020) ---
reg_out <- table1(dep.var.name = "ExcessSupportTrump",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR",
                  reg.data = LongRunData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_C_3D.tex")
write(reg_out, outfile)

rm(LongRunData)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_C_3_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table C.4 ------------------------------------------------------
# The Characteristics of Farmer and Non-Farmer Migrant Households

MigrantsData[, birthOrder := as.numeric(as.character(birthOrder))]
MigrantsData[, numChildren := max(birthOrder, na.rm = T), by = "familyFE"]
MigrantsHH<- MigrantsData[, .SD[1], by = "familyFE", 
                          .SDcols = c("numChildren", "farmerDad", "meanLNIpre",
                                      "ICM", "urban", "f_age")]

MigrantsHH[farmerDad==1, farmerDadW := "Farmers'"]
MigrantsHH[farmerDad==0, farmerDadW := "Non-Farmers'"]

MigrantsHH <- apply_labels(MigrantsHH, 
                           f_age ="Father's Age", 
                           urban ="Urban location", 
                           meanLNIpre ="Avg. Pre-migration LNI score", 
                           numChildren ="Number of Children (0-10)", 
                           ICM = "Intra-Community Marriage")

descTable <- compareGroups(farmerDadW ~ f_age + numChildren + urban + meanLNIpre + ICM, 
                           data = MigrantsHH)

descTableOut <- createTable(descTable, show.p.overall = T)
t <- export2latex(descTableOut, header.labels = c("p.overall" = "\\makecell{Difference\\\\p-value}")) %>% str_split("\n") %>% unlist()

# Format
tt_1 <- regmatches(t[5], gregexpr("[[:digit:]]+", t[5])) %>% unlist() %>% as.numeric() %>% formatC(format="f", digits = 0, big.mark = ",")
tt_0 <- regmatches(t[5], gregexpr("[[:digit:]]+", t[5])) %>% unlist() 
t[5] <- gsub(tt_0[1], tt_1[1], t[5])
t[5] <- gsub(tt_0[2], tt_1[2], t[5])

outfile <- paste0(AppendixTablesPath, "/table_C_4.tex")
write(c("& \\multicolumn{3}{c}{Household:} \\\\",
        "\\\\[-1.8ex]", t[4], "\\\\[-1.8ex]", t[5],
        "\\\\[-1.8ex] & (1) & (2) & (3) \\\\", "\\cline{2-4}", "\\\\[-1.8ex]", 
        "\\\\[-1.8ex]", t[23],
        "\\\\[-1.8ex]", t[24],
        "\\\\[-1.8ex]", t[25],
        "\\\\[-1.8ex]", t[26],
        "\\\\[-1.8ex]", t[27]), outfile)

rm(t, descTable, descTableOut, MigrantsHH)


# Appendix Table C.5 ------------------------------------------------------
# DD and DDD: only rural locations

reg_out <- table3(reg.data = MigrantsData[urban==0])
outfile <- paste0(AppendixTablesPath, "/table_C_5.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Family Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Relative YOB Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_C_5_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table C.6 ------------------------------------------------------
# Selective out-migration: age and children

m_age <- median(OutMigrationData[, age], na.rm = T)
OutMigrationData[, Old := 0]
OutMigrationData[age > m_age, Old := 1]

m_numChild <- median(OutMigrationData[, numChildren], na.rm = T)
OutMigrationData[, manyChild := 0]
OutMigrationData[numChildren > m_numChild, manyChild := 1]


m1 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData)

m2 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData[Old==0])

m3 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData[Old==1])

m4 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData[manyChild==0])

m5 <- felm(formula(paste0("percent ~ SHI25*farmer_y + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData[manyChild==1])

reg.table <- stargazer(m1,m2,m3,m4,m5,
                       float = FALSE,
                       header = FALSE,
                       model.names = FALSE,
                       model.numbers = FALSE,
                       keep = c("SHI25", "farmer_y", "SHI25:farmer_y"),
                       keep.stat = c("n", "rsq"),
                       table.layout = "tas")

reg.table <- str_replace(reg.table, "SHI25:farmer\\\\_y", 
                         "\\\\multirow[t]{2}{*}{\\\\makecell[l]{Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer}}")
reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
reg.table <- str_replace(reg.table, "farmer\\\\_y", "Farmer")
  
outfile <- paste0(AppendixTablesPath, "/table_C_6.tex")
write(c(reg.table[3:10],
        "  & & & & & \\\\ ", "\\hline \\\\[-1.8ex]  \\\\[-1.8ex] ",
        reg.table[12:13]), outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_C_6_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table C.7 ------------------------------------------------------
# Selective out-migration: age and children, interaction

m1 <- felm(formula(paste0("percent ~ SHI25*farmer_y*highLNI + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData)

m2 <- felm(formula(paste0("percent ~ SHI25*farmer_y*highLNI + SHI25*farmer_y*Old + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData)

m3 <- felm(formula(paste0("percent ~ SHI25*farmer_y*highLNI + SHI25*farmer_y*manyChild + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData)

m4 <- felm(formula(paste0("percent ~ SHI25*farmer_y*highLNI + SHI25*farmer_y*Old + SHI25*farmer_y*manyChild + ", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
           data = OutMigrationData)


reg.table <- stargazer(m1,m2,m3,m4,
                       float = FALSE,
                       header = FALSE,
                       model.names = FALSE,
                       model.numbers = FALSE,
                       keep = c("^SHI25$", "^farmer_y$", "^SHI25:farmer_y$", 
                                "^SHI25:farmer_y:highLNI$",
                                "^SHI25:farmer_y:manyChild$",
                                "^SHI25:farmer_y:Old$"),
                       keep.stat = c("n", "rsq"),
                       table.layout = "tas")

reg.table <- str_replace(reg.table, "SHI25:farmer\\\\_y:highLNI", 
                         "\\\\multirow[t]{4}{*}{\\\\makecell[l]{\\\\\\\\ \\\\\\\\ \\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer $\\\\times$ \\\\\\\\ \\\\qquad High Communal \\\\\\\\ \\\\qquad Identification}}")
reg.table <- str_replace(reg.table, "SHI25:farmer\\\\_y:Old", 
                         "\\\\multirow[t]{3}{*}{\\\\makecell[l]{\\\\\\\\ \\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer $\\\\times$ \\\\\\\\ \\\\qquad Age $>$ median}}")
reg.table <- str_replace(reg.table, "SHI25:farmer\\\\_y:manyChild", 
                         "\\\\multirow[t]{3}{*}{\\\\makecell[l]{\\\\\\\\ \\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer $\\\\times$ \\\\\\\\ \\\\qquad Num. Children $>$ median}}")
reg.table <- str_replace(reg.table, "SHI25:farmer\\\\_y", 
                         "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer}}")
reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
reg.table <- str_replace(reg.table, "farmer\\\\_y", "Farmer")

outfile <- paste0(AppendixTablesPath, "/table_C_7.tex")
write(c(reg.table[3:14],
        "  & & & & \\\\ ", 
        "  & & & & \\\\ ", 
        reg.table[15:17], 
        "\\\\", 
        reg.table[18:20], 
        "\\\\[-1.4ex]", 
        "\\hline \\\\[-1.8ex]  \\\\[-1.8ex] ", 
        reg.table[21:22]), outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\")), 
    collapse = "\n") 

outfile <- paste0(AppendixTablesPath, "/table_C_7_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table C.8 ------------------------------------------------------
# Climatic risk 
table_C_8 <- function(dep.var.name1, dep.var.name2){
  
  m1 <- felm(formula(paste0(dep.var.name1, " ~ ClimaticRisk + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  m2 <- felm(formula(paste0(dep.var.name1, " ~ SHI25 + ClimaticRisk + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  m3 <- felm(formula(paste0(dep.var.name1, " ~ SHI25*HighClimaticRisk + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  if(!is.na(dep.var.name2)){
    m4 <- felm(formula(paste0(dep.var.name2, " ~ ClimaticRisk + ", my.SmoothLocCtrls,
                              " | STATE_TERR:YEAR | 0 | grid100m")),
               data = CountyLevelData)
    
    m5 <- felm(formula(paste0(dep.var.name2, " ~ SHI25 + ClimaticRisk + ", my.SmoothLocCtrls,
                              " | STATE_TERR:YEAR | 0 | grid100m")),
               data = CountyLevelData)
    
    m6 <- felm(formula(paste0(dep.var.name2, " ~ SHI25*HighClimaticRisk + ", my.SmoothLocCtrls,
                              " | STATE_TERR:YEAR | 0 | grid100m")),
               data = CountyLevelData)
    
    reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                           float = FALSE,
                           header = FALSE,
                           model.names = FALSE,
                           model.numbers = FALSE,
                           keep = c("SHI25", "ClimaticRisk", 
                                    "HighClimaticRisk"),
                           keep.stat = c("n", "rsq"),
                           table.layout = "tas")
  } else {
    reg.table <- stargazer(m1,m2,m3,
                           float = FALSE,
                           header = FALSE,
                           model.names = FALSE,
                           model.numbers = FALSE,
                           keep = c("SHI25", "ClimaticRisk", 
                                    "HighClimaticRisk"),
                           keep.stat = c("n", "rsq"),
                           table.layout = "tas")
  }
  

  reg.table <- str_replace(reg.table, "SHI25:HighClimaticRisk", 
                           "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad High Climatic Risk}}")
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  reg.table <- str_replace(reg.table, "ClimaticRisk", "Climatic Risk")
  reg.table <- str_replace(reg.table, "HighClimaticRisk", "High Climatic Risk")
  return(c(reg.table[3:13], "\\\\ \\\\[-1.8ex] ", reg.table[c(15:16)]))
}

#--- Panels A and B: LNI and ICM ---
reg_out <- table_C_8(dep.var.name1 = "LNI_county_Native",
                     dep.var.name2 = "ICM_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_8AB.tex")
write(reg_out, outfile)


#--- Panels C and D: TNI and RHI ---
reg_out <- table_C_8(dep.var.name1 = "TNI_Native",
                     dep.var.name2 = "RHI")

outfile <- paste0(AppendixTablesPath, "/table_C_8CD.tex")
write(reg_out, outfile)

#--- Panels E and F: MKP and KPR ---
reg_out <- table_C_8(dep.var.name1 = "MKP_Native",
                     dep.var.name2 = "KPR_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_8EF.tex")
write(reg_out, outfile)


#--- Panel E: SFTI ---
reg_out <- table_C_8(dep.var.name1 = "SFTI_Native",
                     dep.var.name2 = NA)

outfile <- paste0(AppendixTablesPath, "/table_C_8G.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\")), 
    collapse = "\n") 

outfile <- paste0(AppendixTablesPath, "/table_C_8_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table C.9 ------------------------------------------------------
# Trade connectivity: Railroads

table_C_9 <- function(dep.var.name) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(RR)])
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + RR + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(RR)])
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + RR + RR_10m + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(RR)])
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + RR + RR_10m + RR_20m + ", 
                            my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(RR)])
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + RR + RR_10m + RR_20m + RR_30m + ", 
                            my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(RR)])
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + RR + RR_10m + RR_20m + RR_30m + RR_40m + ", 
                            my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(RR)])
  
  m7 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + RR + RR_10m + RR_20m + RR_30m + RR_40m + RR_50m + ", 
                            my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(RR)])
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  
  return(c(reg.table[3:4], "\\\\[-1.8ex] ", reg.table[6:7]))

}


#--- Panel A: LNI ---
reg_out <- table_C_9(dep.var.name = "LNI_county_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_9A.tex")
write(reg_out, outfile)

#--- Panel B: ICM ---
reg_out <- table_C_9(dep.var.name = "ICM_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_9B.tex")
write(reg_out, outfile)

#--- Panel C: TNI ---
reg_out <- table_C_9(dep.var.name = "TNI_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_9C.tex")
write(reg_out, outfile)

#--- Panel D: RHI ---
reg_out <- table_C_9(dep.var.name = "RHI")

outfile <- paste0(AppendixTablesPath, "/table_C_9D.tex")
write(reg_out, outfile)

#--- Panel E: MKP ---
reg_out <- table_C_9(dep.var.name = "MKP_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_9E.tex")
write(reg_out, outfile)

#--- Panel F: KPR ---
reg_out <- table_C_9(dep.var.name = "KPR_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_9F.tex")
write(reg_out, outfile)

#--- Panel G: SFTI ---
reg_out <- table_C_9(dep.var.name = "SFTI_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_9G.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Railroad:", paste(rep(" ", 8), collapse = " & "), "\\\\"),
      paste("\\ \\ \\ \\ In county", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("\\ \\ \\ \\ Within 10 miles", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("\\ \\ \\ \\ Within 20 miles", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("\\ \\ \\ \\ Within 30 miles", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("\\ \\ \\ \\ Within 40 miles", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("\\ \\ \\ \\ Within 50 miles", paste(rep(" ", 8), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n") 

outfile <- paste0(AppendixTablesPath, "/table_C_9_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table C.10 ------------------------------------------------------
# Immigration

table_C_10 <- function(dep.var.name) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + shareImmigrants + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + BirthPlaceDiv + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR + bpl | 0 | grid100m")),
             data = CountyLevelData)
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + shareImmigrants + BirthPlaceDiv + ", 
                            my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR + bpl | 0 | grid100m")),
             data = CountyLevelData)
  
  reg.table <- stargazer(m1,m2,m3,m4,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  
  return(c(reg.table[3:4], "\\\\[-1.8ex] ", reg.table[c(6:7)]))
  
}


#--- Panel A: LNI ---
reg_out <- table_C_10(dep.var.name = "LNI_county_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_10A.tex")
write(reg_out, outfile)


#--- Panel B: ICM ---
reg_out <- table_C_10(dep.var.name = "ICM_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_10B.tex")
write(reg_out, outfile)

#--- Panel C: TNI ---
reg_out <- table_C_10(dep.var.name = "TNI_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_10C.tex")
write(reg_out, outfile)

#--- Panel D: RHI ---
reg_out <- table_C_10(dep.var.name = "RHI")

outfile <- paste0(AppendixTablesPath, "/table_C_10D.tex")
write(reg_out, outfile)

#--- Panel E: MKP ---
reg_out <- table_C_10(dep.var.name = "MKP_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_10E.tex")
write(reg_out, outfile)

#--- Panel E: KPR ---
reg_out <- table_C_10(dep.var.name = "KPR_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_10F.tex")
write(reg_out, outfile)

#--- Panel E: SFTI ---
reg_out <- table_C_10(dep.var.name = "SFTI_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_10G.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- list(c("State $\\times$ Year Fixed Effects", rep("\\checkmark", 4)),
                   c("Geoclimatic Controls", rep("\\checkmark", 4)),
                   c("Smooth Location Controls", rep("\\checkmark", 4)),
                   c("Share Immigrants", "\\checkmark", "", "", "\\checkmark"),
                   c("Birthplace Diversity", "", "\\checkmark", "", "\\checkmark"),
                   c("Dominate Origin Fixed Effects", "", "", "\\checkmark", "\\checkmark")) 
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Share Immigrants & \\checkmark & & & \\checkmark \\\\"),
      paste("Birthplace Diversity & & \\checkmark & & \\checkmark \\\\"),
      paste("Dominate Origin Fixed Effects", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\")), 
    collapse = "\n") 

outfile <- paste0(AppendixTablesPath, "/table_C_10_foot.tex")
write(unlist(footLines), outfile)



# Appendix Table C.11 ------------------------------------------------------
# Race and Slavery 

table_C_11 <- function(dep.var.name1, dep.var.name2){
  
  m1 <- felm(formula(paste0(dep.var.name1, " ~ SHI25 + black_share + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  m2 <- felm(formula(paste0(dep.var.name1, " ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = SlaveData)
  
  m3 <- felm(formula(paste0(dep.var.name1, " ~ SHI25 + slave_share + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = SlaveData)
  
  if(!is.na(dep.var.name2)){
    m4 <- felm(formula(paste0(dep.var.name2, " ~ SHI25 + black_share + ", my.SmoothLocCtrls,
                              " | STATE_TERR:YEAR | 0 | grid100m")),
               data = CountyLevelData)
    
    m5 <- felm(formula(paste0(dep.var.name2, " ~ SHI25 + ", my.SmoothLocCtrls,
                              " | STATE_TERR:YEAR | 0 | grid100m")),
               data = SlaveData)
    
    m6 <- felm(formula(paste0(dep.var.name2, " ~ SHI25 + slave_share +", my.SmoothLocCtrls,
                              " | STATE_TERR:YEAR | 0 | grid100m")),
               data = SlaveData)
    
    reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                           float = FALSE,
                           header = FALSE,
                           model.names = FALSE,
                           model.numbers = FALSE,
                           keep = c("SHI25"),
                           keep.stat = c("n", "rsq"),
                           table.layout = "tas")
  } else {
    reg.table <- stargazer(m1,m2,m3,
                           float = FALSE,
                           header = FALSE,
                           model.names = FALSE,
                           model.numbers = FALSE,
                           keep = c("SHI25"),
                           keep.stat = c("n", "rsq"),
                           table.layout = "tas")
  }
  
  
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  return(c(reg.table[3:4], "\\\\[-1.8ex] ", reg.table[c(6:7)]))
}

FilePath <- paste0(DataPath, "/final/SlaveShareData.csv")
SlaveData <- fread(FilePath, colClasses = my.colClasses)

#--- Panels A and B: LNI and ICM ---
reg_out <- table_C_11(dep.var.name1 = "LNI_county_Native",
                      dep.var.name2 = "ICM_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_11AB.tex")
write(reg_out, outfile)


#--- Panels C and D: TNI and RHI ---
reg_out <- table_C_11(dep.var.name1 = "TNI_Native",
                      dep.var.name2 = "RHI")

outfile <- paste0(AppendixTablesPath, "/table_C_11CD.tex")
write(reg_out, outfile)

#--- Panels E and F: MKP and KPR ---
reg_out <- table_C_11(dep.var.name1 = "MKP_Native",
                      dep.var.name2 = "KPR_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_11EF.tex")
write(reg_out, outfile)

#--- Panel G: SFTI ---
reg_out <- table_C_11(dep.var.name1 = "SFTI_Native",
                      dep.var.name2 = NA)

outfile <- paste0(AppendixTablesPath, "/table_C_11G.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Black Share & \\checkmark & & & \\checkmark & & \\\\"),
      paste("Historical Slave Share & & & \\checkmark & & & \\checkmark \\\\")), 
    collapse = "\n") 

outfile <- paste0(AppendixTablesPath, "/table_C_11_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table C.12 -----------------------------------------------------
# Agricultural inequality 

m1 <- felm(formula(paste0("LNI_county_Native ~ SHI25 + FarmGini +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m2 <- felm(formula(paste0("ICM_Native ~ SHI25 + FarmGini +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m3 <- felm(formula(paste0("TNI_Native ~ SHI25 + FarmGini +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m4 <- felm(formula(paste0("RHI ~ SHI25 + FarmGini +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m5 <- felm(formula(paste0("MKP_Native ~ SHI25 + FarmGini +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m6 <- felm(formula(paste0("KPR_Native ~ SHI25 + FarmGini +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m7 <- felm(formula(paste0("SFTI_Native ~ SHI25 + FarmGini +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                       float = FALSE,
                       header = FALSE,
                       model.names = FALSE,
                       model.numbers = FALSE,
                       keep = c("SHI25"),
                       keep.stat = c("n", "rsq"),
                       table.layout = "tas")

reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")

outfile <- paste0(AppendixTablesPath, "/table_C_12.tex")
write(c(reg.table[3:5], "\\hline \\\\[-1.8ex]  \\\\[-1.8ex] ", reg.table[6:7]), outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Farm Size Gini", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_C_12_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table C.13 -----------------------------------------------------
# Modernization 

table_C_13 <- function(dep.var.name) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(shareUrb2500) & !is.na(shareLitNative) & !is.na(ShareFarmers_Native) & !is.na(ManuEstPC)])

  m2 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + shareUrb2500 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(shareUrb2500) & !is.na(shareLitNative) & !is.na(ShareFarmers_Native) & !is.na(ManuEstPC)])
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + shareLitNative + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(shareUrb2500) & !is.na(shareLitNative) & !is.na(ShareFarmers_Native) & !is.na(ManuEstPC)])
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ShareFarmers_Native + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(shareUrb2500) & !is.na(shareLitNative) & !is.na(ShareFarmers_Native) & !is.na(ManuEstPC)])
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ManuEstPC + ", 
                            my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(shareUrb2500) & !is.na(shareLitNative) & !is.na(ShareFarmers_Native) & !is.na(ManuEstPC)])
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + shareUrb2500 + shareLitNative + ShareFarmers_Native + ManuEstPC + ", 
                            my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData[!is.na(shareUrb2500) & !is.na(shareLitNative) & !is.na(ShareFarmers_Native) & !is.na(ManuEstPC)])
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")

  return(c(reg.table[3:4], "\\\\[-1.8ex] ", reg.table[6:7]))
}

#--- Panel A: LNI ---
reg_out <- table_C_13(dep.var.name = "LNI_county_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_13A.tex")
write(reg_out, outfile)

#--- Panel B: ICM ---
reg_out <- table_C_13(dep.var.name = "ICM_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_13B.tex")
write(reg_out, outfile)

#--- Panel C: TNI ---
reg_out <- table_C_13(dep.var.name = "TNI_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_13C.tex")
write(reg_out, outfile)

#--- Panel D: RHI ---
reg_out <- table_C_13(dep.var.name = "RHI")

outfile <- paste0(AppendixTablesPath, "/table_C_13D.tex")
write(reg_out, outfile)

#--- Panel E: MKP ---
reg_out <- table_C_13(dep.var.name = "MKP_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_13E.tex")
write(reg_out, outfile)

#--- Panel F: KPR ---
reg_out <- table_C_13(dep.var.name = "KPR_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_13F.tex")
write(reg_out, outfile)

#--- Panel G: SFTI ---
reg_out <- table_C_13(dep.var.name = "SFTI_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_13G.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Share Urban & & \\checkmark & & & & \\checkmark \\\\"),
      paste("Share Literate & & & \\checkmark & & & \\checkmark \\\\"),
      paste("Share Farmers & & & & \\checkmark & & \\checkmark \\\\"),
      paste("Manufacturing Est. & & & & & \\checkmark & \\checkmark \\\\")), 
    collapse = "\n")

outfile <- paste0(AppendixTablesPath, "/table_C_13_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table C.14 -----------------------------------------------------
# Frontier experience 

table_C_14 <- function(dep.var.name) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = FrontierStatus)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + frontier100kmL6 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = FrontierStatus)
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR | 0 | grid100m")),
             data = TFE)
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + tye_tfe890_500kNI_100_l6 + ", my.SmoothLocCtrls,
                            " | STATE_TERR | 0 | grid100m")),
             data = TFE)
  
  
  reg.table <- stargazer(m1,m2,m3,m4,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  
  return(c(reg.table[3:4], "\\\\[-1.8ex] ", reg.table[6:7]))
  
}

FilePath <- paste0(DataPath, "/final/CountyLevelDataFrontierStatus.csv")
FrontierStatus <- fread(FilePath, colClasses = my.colClasses)

FilePath <- paste0(DataPath, "/final/CountyLevelData1940TFE.csv")
TFE <- fread(FilePath)

#--- Panel A: LNI ---
reg_out <- table_C_14(dep.var.name = "LNI_county_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_14A.tex")
write(reg_out, outfile)

#--- Panel B: ICM ---
reg_out <- table_C_14(dep.var.name = "ICM_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_14B.tex")
write(reg_out, outfile)

#--- Panel C: TNI ---
reg_out <- table_C_14(dep.var.name = "TNI_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_14C.tex")
write(reg_out, outfile)

#--- Panel D: RHI ---
reg_out <- table_C_14(dep.var.name = "RHI")

outfile <- paste0(AppendixTablesPath, "/table_C_14D.tex")
write(reg_out, outfile)

#--- Panel E: MKP ---
reg_out <- table_C_14(dep.var.name = "MKP_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_14E.tex")
write(reg_out, outfile)

#--- Panel F: KPR ---
reg_out <- table_C_14(dep.var.name = "KPR_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_14F.tex")
write(reg_out, outfile)

#--- Panel G: SFTI ---
reg_out <- table_C_14(dep.var.name = "SFTI_Native")

outfile <- paste0(AppendixTablesPath, "/table_C_14G.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects & \\checkmark & \\checkmark & & \\\\"),
      paste("State Fixed Effects & & & \\checkmark & \\checkmark \\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Frontier Status & & \\checkmark & & \\\\"),
      paste("Total Frontier Experience & & & & \\checkmark \\\\")), 
    collapse = "\n")

outfile <- paste0(AppendixTablesPath, "/table_C_14_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.1 ------------------------------------------------------
# Drop SHI outliers

table_D_1 <- function(data.restriction) {
  
  m1 <- felm(formula(paste0("LNI_county_Native ~ SHI25 +", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData %>% filter(eval(parse(text = data.restriction))))
  
  m2 <- felm(formula(paste0("ICM_Native ~ SHI25 +", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData %>% filter(eval(parse(text = data.restriction))))
  
  m3 <- felm(formula(paste0("TNI_Native ~ SHI25 +", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData %>% filter(eval(parse(text = data.restriction))))
  
  m4 <- felm(formula(paste0("RHI ~ SHI25 +", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData %>% filter(eval(parse(text = data.restriction))))
  
  m5 <- felm(formula(paste0("MKP_Native ~ SHI25 +", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData %>% filter(eval(parse(text = data.restriction))))
  
  m6 <- felm(formula(paste0("KPR_Native ~ SHI25 +", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData %>% filter(eval(parse(text = data.restriction))))
  
  m7 <- felm(formula(paste0("SFTI_Native ~ SHI25 +", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData %>% filter(eval(parse(text = data.restriction))))
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = c("SHI25"),
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  
  return(reg.table[3:7])
  
}

#--- Panel A: 0.25 < SHI < 0.75 ---
reg_out <- table_D_1(data.restriction = "SHI25 > .25 & SHI25 < .75")

outfile <- paste0(AppendixTablesPath, "/table_D_1A.tex")
write(reg_out, outfile)


#--- Panel B: p(5) < SHI < p(95) ---
q05 <- quantile(CountyLevelData$SHI25, 0.05, na.rm=T)
q95 <- quantile(CountyLevelData$SHI25, 0.95, na.rm=T)

reg_out <- table_D_1(data.restriction = "SHI25 > q05 & SHI25 < q95")

outfile <- paste0(AppendixTablesPath, "/table_D_1B.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_1_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.2 ------------------------------------------------------
# Drop state FEs

#--- Panel A: LNI ---
reg_out <- table1(dep.var.name = "LNI_county_Native",
                  indep.var.name = "SHI25",
                  FEs = "YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_2A.tex")
write(reg_out, outfile)

#--- Panel B: ICM ---
reg_out <- table1(dep.var.name = "ICM_Native",
                  indep.var.name = "SHI25",
                  FEs = "YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_2B.tex")
write(reg_out, outfile)

#--- Panel C: TNI ---
reg_out <- table1(dep.var.name = "TNI_Native",
                  indep.var.name = "SHI25",
                  FEs = "YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_2C.tex")
write(reg_out, outfile)

#--- Panel D: RHI ---
reg_out <- table1(dep.var.name = "RHI",
                  indep.var.name = "SHI25",
                  FEs = "YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_2D.tex")
write(reg_out, outfile)

#--- Panel E: MKP ---
reg_out <- table1(dep.var.name = "MKP_Native",
                  indep.var.name = "SHI25",
                  FEs = "YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_2E.tex")
write(reg_out, outfile)

#--- Panel F: KPR ---
reg_out <- table1(dep.var.name = "KPR_Native",
                  indep.var.name = "SHI25",
                  FEs = "YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_2F.tex")
write(reg_out, outfile)

#--- Panel G: SFTI ---
reg_out <- table1(dep.var.name = "SFTI_Native",
                  indep.var.name = "SHI25",
                  FEs = "YEAR",
                  reg.data = CountyLevelData,
                  clusterVar = "grid100m",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_2G.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_2_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.3 ------------------------------------------------------
# Different Definitions of LNI

m1 <- felm(formula(paste0("LNI_county_Native ~ SHI25 +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m2 <- felm(formula(paste0("LNI_county_NativeNYSIIS ~ SHI25 +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m3 <- felm(formula(paste0("LNI_county_Native100 ~ SHI25 +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m4 <- felm(formula(paste0("LNI_county_Native1000 ~ SHI25 +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

m5 <- felm(formula(paste0("LNI_county_Native ~ SHI25 +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData[childNNative >= 100])

m6 <- felm(formula(paste0("LNI_state_Native ~ SHI25 +", my.SmoothLocCtrls,
                          " | STATE_TERR:YEAR | 0 | grid100m")),
           data = CountyLevelData)

reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                       float = FALSE,
                       header = FALSE,
                       model.names = FALSE,
                       model.numbers = FALSE,
                       keep = c("SHI25"),
                       keep.stat = c("n", "rsq"),
                       add.lines = c(list(c("Dependent Variable Mean", 
                                            mean(m1$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            mean(m2$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            mean(m3$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            mean(m4$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            mean(m5$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            mean(m6$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ",")),
                                          c("Dependent Variable SD", 
                                            sd(m1$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            sd(m2$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            sd(m3$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            sd(m4$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            sd(m5$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ","),
                                            sd(m6$response) %>% 
                                              formatC(format="f", digits = 2, big.mark = ",")))),
                       table.layout = "tas")

reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")

outfile <- paste0(AppendixTablesPath, "/table_D_3.tex")
write(c(reg.table[3:4], "& & & & \\\\ \\hline \\\\[-1.8ex]  \\\\[-1.8ex] ", reg.table[6:9]), outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_3_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.4 ------------------------------------------------------
# Different census samples (immigrants and non-Whites)

table_D_4 <- function(dep.var.name) {
  
  m1 <- felm(formula(paste0(dep.var.name, "_Native ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ SHI25  + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  m3 <- felm(formula(paste0(dep.var.name, "_NativeAllRace ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  m4 <- felm(formula(paste0(dep.var.name, "_AllRace ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = CountyLevelData)
  
  
  reg.table <- stargazer(m1,m2,m3,m4,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         add.lines = c(list(c("Dependent Variable Mean", 
                                              mean(m1$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m2$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m3$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              mean(m4$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ",")),
                                            c("Dependent Variable SD", 
                                              sd(m1$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m2$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m3$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ","),
                                              sd(m4$response) %>% 
                                                formatC(format="f", digits = 2, big.mark = ",")))),table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")
  
  return(c(reg.table[3:4], "\\\\[-1.8ex]", reg.table[6:9]))
}

#--- Panel A: LNI ---
reg_out <- table_D_4(dep.var.name = "LNI_county")

outfile <- paste0(AppendixTablesPath, "/table_D_4A.tex")
write(reg_out, outfile)

#--- Panel B: ICM ---
reg_out <- table_D_4(dep.var.name = "ICM")

outfile <- paste0(AppendixTablesPath, "/table_D_4B.tex")
write(reg_out, outfile)

#--- Panel C: TNI ---
reg_out <- table_D_4(dep.var.name = "TNI")

outfile <- paste0(AppendixTablesPath, "/table_D_4C.tex")
write(reg_out, outfile)

#--- Panel D: MKP ---
reg_out <- table_D_4(dep.var.name = "MKP")

outfile <- paste0(AppendixTablesPath, "/table_D_4D.tex")
write(reg_out, outfile)

#--- Panel E: KPR ---
reg_out <- table_D_4(dep.var.name = "KPR")

outfile <- paste0(AppendixTablesPath, "/table_D_4E.tex")
write(reg_out, outfile)

#--- Panel G: SFTI ---
reg_out <- table_D_4(dep.var.name = "SFTI")

outfile <- paste0(AppendixTablesPath, "/table_D_4F.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\")), 
    collapse = "\n") 

outfile <- paste0(AppendixTablesPath, "/table_D_4_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.5 ------------------------------------------------------
# Different distances for SHI calculation 

table_D_5 <- function(dep.var.name, reg.data) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ SHI25 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = reg.data)
  rownames(m1$coefficients)[1] <- "Soil Heterogeneity"
  rownames(m1$beta)[1] <- "Soil Heterogeneity"  
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ SHI10 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = reg.data)
  rownames(m2$coefficients)[1] <- "Soil Heterogeneity"
  rownames(m2$beta)[1] <- "Soil Heterogeneity"  
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ SHI15 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = reg.data)
  rownames(m3$coefficients)[1] <- "Soil Heterogeneity"
  rownames(m3$beta)[1] <- "Soil Heterogeneity"  
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ SHI20 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = reg.data)
  rownames(m4$coefficients)[1] <- "Soil Heterogeneity"
  rownames(m4$beta)[1] <- "Soil Heterogeneity"  
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ SHI30 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = reg.data)
  rownames(m5$coefficients)[1] <- "Soil Heterogeneity"
  rownames(m5$beta)[1] <- "Soil Heterogeneity"  
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ SHI35 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = reg.data)
  rownames(m6$coefficients)[1] <- "Soil Heterogeneity"
  rownames(m6$beta)[1] <- "Soil Heterogeneity"  
  
  m7 <- felm(formula(paste0(dep.var.name, " ~ SHI40 + ", my.SmoothLocCtrls,
                            " | STATE_TERR:YEAR | 0 | grid100m")),
             data = reg.data)
  rownames(m7$coefficients)[1] <- "Soil Heterogeneity"  
  rownames(m7$beta)[1] <- "Soil Heterogeneity"  
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "Soil Heterogeneity",
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  return(c(reg.table[3:4], "\\\\[-1.8ex]", reg.table[6:7]))
  
}

#--- Panel A: LNI ---
reg_out <- table_D_5(dep.var.name = "LNI_county_Native", reg.data = CountyLevelData)

outfile <- paste0(AppendixTablesPath, "/table_D_5A.tex")
write(reg_out, outfile)

#--- Panel B: ICM ---
reg_out <- table_D_5(dep.var.name = "ICM_Native", reg.data = CountyLevelData)

outfile <- paste0(AppendixTablesPath, "/table_D_5B.tex")
write(reg_out, outfile)

#--- Panel C: TNI ---
reg_out <- table_D_5(dep.var.name = "TNI_Native", reg.data = CountyLevelData)

outfile <- paste0(AppendixTablesPath, "/table_D_5C.tex")
write(reg_out, outfile)

#--- Panel D: RHI ---
reg_out <- table_D_5(dep.var.name = "RHI", reg.data = CountyLevelData)

outfile <- paste0(AppendixTablesPath, "/table_D_5D.tex")
write(reg_out, outfile)

#--- Panel E: MKP ---
reg_out <- table_D_5(dep.var.name = "MKP_Native", reg.data = CountyLevelData)

outfile <- paste0(AppendixTablesPath, "/table_D_5E.tex")
write(reg_out, outfile)

#--- Panel F: KPR ---
reg_out <- table_D_5(dep.var.name = "KPR_Native", reg.data = CountyLevelData)

outfile <- paste0(AppendixTablesPath, "/table_D_5F.tex")
write(reg_out, outfile)

#--- Panel G: SFTI ---
reg_out <- table_D_5(dep.var.name = "SFTI_Native", reg.data = CountyLevelData)

outfile <- paste0(AppendixTablesPath, "/table_D_5G.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_5_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table D.6 ------------------------------------------------------
# Selection of migrants - no origin FEs

table_D_6 <- function(dep.var.name, indep.var.name, FEs, ind.ctrls, reg.data) {
  
  #--- Estimate the seven specifications ---
  m1 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, 
                            " | ", FEs, " | 0 | GISJOIN")), 
             data = reg.data)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.GeoCtrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.AgSuitabilityCtrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.HigherOrderCtrls,
                            " | ", FEs, " | 0 | GISJOIN")),
             data = reg.data)
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.HigherOrderCtrls,
                            " | ", FEs, " + ", ind.ctrls, " | 0 | GISJOIN")),
             data = reg.data)
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "SHI25", "Soil Heterogeneity")

  return(c(reg.table[3:5], "\\\\[-1.8ex] ", reg.table[6:7]))
}

#--- Panel A: Farming Experience ---

reg_out <- table_D_6(dep.var.name = "farmer_z",
                  indep.var.name = "SHI25",
                  FEs = "STATE_TERR_y:YEAR",
                  ind.ctrls = "age_y",
                  reg.data = SelectionFarmingData)

outfile <- paste0(AppendixTablesPath, "/table_D_6A.tex")
write(reg_out, outfile)


#--- Panel B: Children's LNI ---

reg_out <- table_D_6(dep.var.name = "LNI_state_Native",
                     indep.var.name = "SHI25",
                     FEs = "statefips:YEAR",
                     ind.ctrls = "sex:birthYear + birthOrder",
                     reg.data = MigrantsData[Post==0])

outfile <- paste0(AppendixTablesPath, "/table_D_6B.tex")
write(reg_out, outfile)

#--- Panel C: ICM ---
reg_out <- table_D_6(dep.var.name = "ICM",
                  indep.var.name = "SHI25",
                  FEs = "statefips:YEAR",
                  ind.ctrls = "sex:birthYear + f_age",
                  reg.data = MigrantsData[, .SD[1], by = 'familyFE'])

outfile <- paste0(AppendixTablesPath, "/table_D_6C.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"), 
      paste("Individual Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_6_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table D.7 ------------------------------------------------------
# DD and DDD, robustness - individual controls 


reg_out <- table3(reg.data = MigrantsData, ind.ctrls = "sex + birthOrder + cohortFE5Y")
outfile <- paste0(AppendixTablesPath, "/table_D_7.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Family Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Relative YOB Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Individual Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_7_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.8 ------------------------------------------------------
# DD and DDD, robustness - samples (immigrants and non-Whites)

table_D_8 <- function(reg.data, dep.var.name) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ Post:SHI25 | familyFE + Rtime | 0 | GISJOIN")), 
             data = reg.data)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ Post:SHI25 | familyFE + Rtime | 0 | GISJOIN")), 
             data = reg.data[farmerDad==1])
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ Post:SHI25 | familyFE + Rtime | 0 | GISJOIN")), 
             data = reg.data[farmerDad==0])
  
  reg.table <- stargazer(m1,m2,m3,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = c("Post:SHI25"),
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  reg.table <- str_replace(reg.table, "Post:SHI25", 
                           "\\\\multirow[t]{2}{*}{\\\\makecell[l]{Post Migration  $\\\\times$ \\\\\\\\ \\\\quad Soil Heterogeneity}}")
  
  return(reg.table[3:5])
}
  

#--- Panel A: Include foreign-born parents ---
MyFilePath <- paste0(DataPath, "/final/migrants/Migrants_2.csv")
MigrantsData <- fread(MyFilePath)

reg_out <- table_D_8(reg.data = MigrantsData, dep.var.name = "LNI_state")
outfile <- paste0(AppendixTablesPath, "/table_D_8A.tex")
write(reg_out, outfile)

#--- Panel B: Include all races ---
MyFilePath <- paste0(DataPath, "/final/migrants/MigrantsNativeAllRace_2.csv")
MigrantsData <- fread(MyFilePath)

reg_out <- table_D_8(reg.data = MigrantsData, dep.var.name = "LNI_state_NativeAllRace")
outfile <- paste0(AppendixTablesPath, "/table_D_8B.tex")
write(reg_out, outfile)

#--- Panel C: Include multiple moves ---
MyFilePath <- paste0(DataPath, "/final/migrants/MigrantsNative.csv")
MigrantsData <- fread(MyFilePath)

reg_out <- table_D_8(reg.data = MigrantsData, dep.var.name = "LNI_state_Native")
outfile <- paste0(AppendixTablesPath, "/table_D_8C.tex")
write(reg_out, outfile)

rm(MigrantsData)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Family Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Relative YOB Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_8_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.9 ------------------------------------------------------
# The Long-Run Impact on Communal Ties: Geographical Mobility
# Robustness 1: Controls

#--- Panel A: All Migrants' Children ---
reg_out <- table1(dep.var.name = "percent",
                  indep.var.name = "SHI25",
                  FEs = "bpl:statefip_y:YEAR:Rtime:highLNI:ICM",
                  reg.data = MigrantsChildren,
                  clusterVar = "GISJOIN",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_9A.tex")
write(reg_out, outfile)

#--- Panel B: Farmers' Children ---
reg_out <- table1(dep.var.name = "percent",
                  indep.var.name = "SHI25",
                  FEs = "bpl:statefip_y:YEAR:Rtime:highLNI:ICM",
                  reg.data = MigrantsChildren[farmerDad==1],
                  clusterVar = "GISJOIN",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_9B.tex")
write(reg_out, outfile)

#--- Panel C: All Migrants' Children, interaction ---
reg_out <- table1(dep.var.name = "percent",
                  indep.var.name = "SHI25*farmerDad",
                  FEs = "bpl:statefip_y:YEAR:Rtime:highLNI:ICM",
                  reg.data = MigrantsChildren,
                  clusterVar = "GISJOIN",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25:farmerDad", 
                       "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmers$'$ Household}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7], "\\\\[-1.8ex] ", reg_out[9:10])

outfile <- paste0(AppendixTablesPath, "/table_D_9C.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Cell Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_9_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table D.10 ------------------------------------------------------
# The Long-Run Impact on Communal Ties: Geographical Mobility
# Robustness 2: Cells

table_D_10 <- function(dep.var.name, indep.var.name, reg.data) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | 1 | 0 | GISJOIN")),
             data = reg.data)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | statefip_y | 0 | GISJOIN")),
             data = reg.data)
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | YEAR:statefip_y | 0 | GISJOIN")),
             data = reg.data)
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | bpl:statefip_y:YEAR | 0 | GISJOIN")),
             data = reg.data)
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | bpl:statefip_y:YEAR:Rtime | 0 | GISJOIN")),
             data = reg.data)
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | bpl:statefip_y:YEAR:Rtime:highLNI:ICM | 0 | GISJOIN")),
             data = reg.data)
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = "SHI25",
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  return(reg.table)
}

#--- Panel A: All Migrants' Children ---
reg_out <- table_D_10(dep.var.name = "percent",
                      indep.var.name = "SHI25",
                      reg.data = MigrantsChildren)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

outfile <- paste0(AppendixTablesPath, "/table_D_10A.tex")
write(c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7]), outfile)


#--- Panel B: Farmers' Children ---
reg_out <- table_D_10(dep.var.name = "percent",
                      indep.var.name = "SHI25",
                      reg.data = MigrantsChildren[farmerDad==1])

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

outfile <- paste0(AppendixTablesPath, "/table_D_10B.tex")
write(c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7]), outfile)


#--- Panel C: All Migrants' Children, interaction ---
reg_out <- table_D_10(dep.var.name = "percent",
                      indep.var.name = "SHI25*farmerDad",
                      reg.data = MigrantsChildren)

reg_out <- str_replace(reg_out, "SHI25:farmerDad", 
                         "\\\\multirow[t]{2}{*}{\\\\makecell[l]{Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmers$'$ Household}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7], "\\\\[-1.8ex] ", reg_out[9:10])

outfile <- paste0(AppendixTablesPath, "/table_D_10C.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Destination State", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("\\ $\\times$ Census year", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("\\ $\\times$ Birth State", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("\\ $\\times$ Relative-year-of-birth", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("\\ $\\times$ High Parental Communalism", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_10_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.11 ------------------------------------------------------
# The Long-Run Impact on Communal Ties: ICM
# robustness 1: Controls

#--- Panel A: All Migrants' Children ---
reg_out <- table1(dep.var.name = "ICM_d",
                  indep.var.name = "SHI25",
                  FEs = "bpl:statefip_y:YEAR:Rtime:highLNI:ICM",
                  reg.data = MigrantsChildren,
                  clusterVar = "GISJOIN",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_11A.tex")
write(reg_out, outfile)

#--- Panel B: Farmers' Children ---
reg_out <- table1(dep.var.name = "ICM_d",
                  indep.var.name = "SHI25",
                  FEs = "bpl:statefip_y:YEAR:Rtime:highLNI:ICM",
                  reg.data = MigrantsChildren[farmerDad==1],
                  clusterVar = "GISJOIN",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7])

outfile <- paste0(AppendixTablesPath, "/table_D_11B.tex")
write(reg_out, outfile)


#--- Panel C: All Migrants' Children, interaction ---
reg_out <- table1(dep.var.name = "ICM_d",
                  indep.var.name = "SHI25*farmerDad",
                  FEs = "bpl:statefip_y:YEAR:Rtime:highLNI:ICM",
                  reg.data = MigrantsChildren,
                  clusterVar = "GISJOIN",
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25:farmerDad", 
                       "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmers$'$ Household}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7], "\\\\[-1.8ex] ", reg_out[9:10])

outfile <- paste0(AppendixTablesPath, "/table_D_11C.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Cell Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_11_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.12 ------------------------------------------------------
# The Long-Run Impact on Communal Ties: ICM
# Robustness 2: Cells

#--- Panel A: All Migrants' Children ---
reg_out <- table_D_10(dep.var.name = "ICM_d",
                      indep.var.name = "SHI25",
                      reg.data = MigrantsChildren)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

outfile <- paste0(AppendixTablesPath, "/table_D_12A.tex")
write(c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7]), outfile)


#--- Panel B: Farmers' Children ---
reg_out <- table_D_10(dep.var.name = "ICM_d",
                      indep.var.name = "SHI25",
                      reg.data = MigrantsChildren[farmerDad==1])

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

outfile <- paste0(AppendixTablesPath, "/table_D_12B.tex")
write(c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7]), outfile)


#--- Panel C: All Migrants' Children, interaction ---
reg_out <- table_D_10(dep.var.name = "ICM_d",
                      indep.var.name = "SHI25*farmerDad",
                      reg.data = MigrantsChildren)

reg_out <- str_replace(reg_out, "SHI25:farmerDad", 
                       "\\\\multirow[t]{2}{*}{\\\\makecell[l]{Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmers$'$ Household}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")

outfile <- paste0(AppendixTablesPath, "/table_D_12C.tex")
reg_out <- c(reg_out[3:4], "\\\\[-1.8ex] ", reg_out[6:7], "\\\\[-1.8ex] ", reg_out[9:10])
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Destination State", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("\\ $\\times$ Census year", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("\\ $\\times$ Birth State", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("\\ $\\times$ Relative-year-of-birth", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("\\ $\\times$ High Parental Communalism", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_12_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.13 ------------------------------------------------------
# Robustness of learning results, different measures

#--- Learning potential: 3 crops ---
reg_out <- table5(dep.var.name = "g_shareFertilizer",
                  FEs = "STATE_TERR:YEAR",
                  learning.var = "gains3", 
                  high.learning.var = "high.gains3",
                  reg.ctrls = my.SmoothLocCtrls,
                  reg.data = FertilizerData)

outfile <- paste0(AppendixTablesPath, "/table_D_13A.tex")
write(reg_out, outfile)


#--- Growth in Wheat Production ---
FilePath <- paste0(DataPath, "/final/WheatProductionData.csv")
WheatProductionData <- fread(FilePath, colClasses = my.colClasses)

reg_out <- table5(dep.var.name = "g_wheatProduction",
                  FEs = "STATE_TERR:YEAR",
                  learning.var = "WheatDiff", 
                  high.learning.var = "high.gains",
                  reg.ctrls = my.SmoothLocCtrls,
                  reg.data = WheatProductionData)

outfile <- paste0(AppendixTablesPath, "/table_D_13B.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_13_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table D.14 -----------------------------------------------------------------
# Robustness of growth in fertilizers use, controls


table_D_14 <- function(dep.var.name, indep.var.name, FEs, reg.data, clusterVar,
                       OsterDelta = FALSE) {
  
  #--- Estimate the six specifications ---
  m1 <- felm(formula(paste0("g_", dep.var.name, " ~ ", indep.var.name, " | 1 | 0 | ", clusterVar)), 
             data = reg.data)
  
  m2 <- felm(formula(paste0("g_", dep.var.name, " ~ ", indep.var.name, " | ", FEs, " | 0 | ", clusterVar)), 
             data = reg.data)
  
  m3 <- felm(formula(paste0("g_", dep.var.name, " ~ ", indep.var.name, " + ", my.GeoCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  m4 <- felm(formula(paste0("g_", dep.var.name, " ~ ", indep.var.name, " + ", my.SmoothLocCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  m5 <- felm(formula(paste0("g_", dep.var.name, " ~ ", indep.var.name, " + ", my.AgSuitabilityCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  m6 <- felm(formula(paste0("g_", dep.var.name, " ~ ", indep.var.name, " + ", my.HigherOrderCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  m7 <- felm(formula(paste0("g_", dep.var.name, " ~ ", indep.var.name, " + l_", dep.var.name, " + ", my.HigherOrderCtrls,
                            " | ", FEs, " | 0 | ", clusterVar)),
             data = reg.data)
  
  if(OsterDelta) {
    reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                           float = FALSE,
                           header = FALSE,
                           model.names = FALSE,
                           model.numbers = FALSE,
                           keep = c("SHI25", paste0("l_", dep.var.name)),
                           keep.stat = c("n", "rsq"),
                           add.lines = list(c("Oster $\\delta$ for $\\beta=0$", "", 
                                              o_delta(y = paste0("g_", dep.var.name),
                                                      x = "SHI25",
                                                      con = paste0("+ factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m7)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2), 
                                              o_delta(y = paste0("g_", dep.var.name),
                                                      x = "SHI25",
                                                      con = paste0("+ ", my.GeoCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m7)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2), 
                                              o_delta(y = paste0("g_", dep.var.name),
                                                      x = "SHI25",
                                                      con = paste0("+ ", my.SmoothLocCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m7)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2), 
                                              o_delta(y = paste0("g_", dep.var.name),
                                                      x = "SHI25",
                                                      con = paste0("+ ", my.AgSuitabilityCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m7)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2), 
                                              o_delta(y = paste0("g_", dep.var.name),
                                                      x = "SHI25",
                                                      con = paste0("+ ", my.HigherOrderCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m7)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2),
                                              o_delta(y = paste0("g_", dep.var.name),
                                                      x = "SHI25",
                                                      con = paste0("+ l_", dep.var.name, " + ", my.HigherOrderCtrls, " + factor(", FEs, ")"),
                                                      beta = 0,
                                                      R2max = min(1, 1.3*summary(m7)$r2),
                                                      type = "lm",
                                                      data = reg.data)[1,2] %>%
                                                as.numeric() %>%
                                                formatC(format="f", digits = 2))),
                           table.layout = "tas")
    
    return(reg.table)
  } else{
    reg.table <- stargazer(m1,m2,m3,m4,m5,m6,m7,
                           float = FALSE,
                           header = FALSE,
                           model.names = FALSE,
                           model.numbers = FALSE,
                           keep = c("SHI25", paste0("l_", dep.var.name)),
                           keep.stat = c("n", "rsq"),
                           table.layout = "tas")
    
    return(reg.table)
  }
}

#--- Panel A: All counties ---
reg_out <- table_D_14(dep.var.name = "shareFertilizer",
                      indep.var.name = "SHI25",
                      FEs = "STATE_TERR:YEAR",
                      reg.data = FertilizerData,
                      clusterVar = "grid100m",
                      OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "l\\\\_shareFertilizer", "Share Using Fertilizer$_{t-1}$")

reg_out <- c(reg_out[3:8], "\\\\[-1.8ex] ", reg_out[9:11])

outfile <- paste0(AppendixTablesPath, "/table_D_14A.tex")
write(reg_out, outfile)


#--- Panel B: High Learning Potential Counties ---
reg_out <- table_D_14(dep.var.name = "shareFertilizer",
                      indep.var.name = "SHI25",
                      FEs = "STATE_TERR:YEAR",
                      reg.data = FertilizerData[high.gains5 == 1],
                      clusterVar = "grid100m",
                      OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "l\\\\_shareFertilizer", "Share Using Fertilizer$_{t-1}$")

reg_out <- c(reg_out[3:8], "\\\\[-1.8ex] ", reg_out[9:10])

outfile <- paste0(AppendixTablesPath, "/table_D_14B.tex")
write(reg_out, outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_14_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table D.15 -----------------------------------------------------------------
# Robustness of growth in wheat share in acres, controls

#--- Panel A: All counties ---
reg_out <- table_D_14(dep.var.name = "shareWheat",
                      indep.var.name = "SHI25",
                      FEs = "STATE_TERR:YEAR",
                      reg.data = WheatShareData,
                      clusterVar = "grid100m",
                      OsterDelta = TRUE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "l\\\\_shareWheat", "Wheat Share$_{t-1}$")

reg_out <- c(reg_out[3:8], "\\\\[-1.8ex] ", reg_out[9:11])

outfile <- paste0(AppendixTablesPath, "/table_D_15A.tex")
write(reg_out, outfile)


#--- Panel B: High Learning Potential Counties ---
reg_out <- table_D_14(dep.var.name = "shareWheat",
                      indep.var.name = "SHI25",
                      FEs = "STATE_TERR:YEAR",
                      reg.data = WheatShareData[high.gains == 1],
                      clusterVar = "grid100m",
                      OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "l\\\\_shareWheat", "Wheat Share$_{t-1}$")

reg_out <- c(reg_out[3:8], "\\\\[-1.8ex] ", reg_out[9:10])

outfile <- paste0(AppendixTablesPath, "/table_D_15B.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_18_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table D.16 -----------------------------------------------------------------
# Robustness of interaction with high learning potential, controls

reg_out <- table1(dep.var.name = "g_shareFertilizer",
                  indep.var.name = "SHI25*high.gains5",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = FertilizerData,
                  clusterVar = "grid100m",
                  keep.vars = c("SHI25", "high.gains5"),
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25:high.gains5", 
                       "\\\\multirow[t]{2}{*}{\\\\makecell[l]{Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad High Learning Potential}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "high.gains5", "High Learning Potential")
reg_out <- c(reg_out[3:11], "\\\\[-1.8ex] ", reg_out[12:13])

outfile <- paste0(AppendixTablesPath, "/table_D_16A.tex")
write(reg_out, outfile)



reg_out <- table1(dep.var.name = "g_shareWheat",
                  indep.var.name = "SHI25*high.gains",
                  FEs = "STATE_TERR:YEAR",
                  reg.data = WheatShareData,
                  clusterVar = "grid100m",
                  keep.vars = c("SHI25", "high.gains"),
                  OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25:high.gains", 
                       "\\\\multirow[t]{2}{*}{\\\\makecell[l]{Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad High Learning Potential}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "high.gains", "High Learning Potential")
reg_out <- c(reg_out[3:11], "\\\\[-1.8ex] ", reg_out[12:13])

outfile <- paste0(AppendixTablesPath, "/table_D_16B.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_16_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table D.17 ------------------------------------------------------
# Drop state FEs

#--- Panel A: Fertilizer ---
reg_out <- table_D_14(dep.var.name = "shareFertilizer",
                      indep.var.name = "SHI25",
                      FEs = "YEAR",
                      reg.data = FertilizerData,
                      clusterVar = "grid100m",
                      OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "l\\\\_shareFertilizer", "Share Using Fertilizer$_{t-1}$")

reg_out <- c(reg_out[3:8], "\\\\[-1.8ex] ", reg_out[9:10])

outfile <- paste0(AppendixTablesPath, "/table_D_17A.tex")
write(reg_out, outfile)



#--- Panel B: Wheat Share ---
reg_out <- table_D_14(dep.var.name = "shareWheat",
                      indep.var.name = "SHI25",
                      FEs = "YEAR",
                      reg.data = WheatShareData,
                      clusterVar = "grid100m",
                      OsterDelta = FALSE)

reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "l\\\\_shareWheat", "Wheat Share$_{t-1}$")

reg_out <- c(reg_out[3:8], "\\\\[-1.8ex] ", reg_out[9:10])

outfile <- paste0(AppendixTablesPath, "/table_D_17B.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 6), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_17_foot.tex")
write(unlist(footLines), outfile)

# Appendix Table D.18 -----------------------------------------------------------------
# Robustness of learning results to different distances for SHI calculation

#--- Panel A: Fertilizer ---
reg_out <- table_D_5(dep.var.name = "g_shareFertilizer", reg.data = FertilizerData)

outfile <- paste0(AppendixTablesPath, "/table_D_18A.tex")
write(reg_out, outfile)


#--- Panel B: Wheat Share ---
reg_out <- table_D_5(dep.var.name = "g_shareWheat", reg.data = WheatShareData)

outfile <- paste0(AppendixTablesPath, "/table_D_18B.tex")
write(reg_out, outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 2), collapse = " & "), paste(rep("\\checkmark", 7), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_18_foot.tex")
write(unlist(footLines), outfile)



# Appendix Table D.19 -----------------------------------------------------
# Robustness of selective out-migration 1. High LNI, ctrls

table_D_19 <- function(dep.var.name, indep.var.name, reg.data) {
  
  m1 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, 
                            " | 1 | 0 | GISJOIN_y")),
             data = reg.data)
  
  m2 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, 
                            " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
             data = reg.data)
  
  m3 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.GeoCtrls,
                            " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
             data = reg.data)
  
  m4 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.AgSuitabilityCtrls,
                            " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
             data = reg.data)
  
  m5 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.HigherOrderCtrls,
                            " | STATE_TERR:YEAR | 0 | GISJOIN_y")),
             data = reg.data)
  
  m6 <- felm(formula(paste0(dep.var.name, " ~ ", indep.var.name, " + ", my.HigherOrderCtrls,
                            " | STATE_TERR:YEAR + age + bpl | 0 | GISJOIN_y")),
             data = reg.data)
  
  reg.table <- stargazer(m1,m2,m3,m4,m5,m6,
                         float = FALSE,
                         header = FALSE,
                         model.names = FALSE,
                         model.numbers = FALSE,
                         keep = c("SHI25", "farmer_y"),
                         keep.stat = c("n", "rsq"),
                         table.layout = "tas")
  
  return(reg.table)
}

reg_out <- table_D_19(dep.var.name = "percent",
                      indep.var.name = "SHI25*farmer_y",
                      reg.data = OutMigrationData[highLNI==1])

reg_out <- str_replace(reg_out, "SHI25:farmer\\\\_y", 
                       "\\\\multirow[t]{2}{*}{\\\\makecell[l]{Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "farmer\\\\_y", "Farmer")

outfile <- paste0(AppendixTablesPath, "/table_D_19.tex")
write(c(reg_out[3:11], "\\hline \\\\[-1.8ex]  \\\\[-1.8ex] ", reg_out[12:13]), outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"), 
      paste("Individual Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n")  

outfile <- paste0(AppendixTablesPath, "/table_D_19_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.20----------------------------------------------------
# Robustness of selective out-migration 1. ICM, ctrls

reg_out <- table_D_19(dep.var.name = "percent",
                      indep.var.name = "SHI25*farmer_y",
                      reg.data = OutMigrationData[ICM==1])

reg_out <- str_replace(reg_out, "SHI25:farmer\\\\_y", 
                       "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "farmer\\\\_y", "Farmer")

outfile <- paste0(AppendixTablesPath, "/table_D_20.tex")
write(c(reg_out[3:11], "\\hline \\\\[-1.8ex]  \\\\[-1.8ex] ", reg_out[12:13]), outfile)

#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"), 
      paste("Individual Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n") 

outfile <- paste0(AppendixTablesPath, "/table_D_20_foot.tex")
write(unlist(footLines), outfile)


# Appendix Table D.21 -----------------------------------------------------
# Robustness of selective out-migration 2. Interaction with high LNI, ctrls


reg_out <- table_D_19(dep.var.name = "percent",
                      indep.var.name = "SHI25*farmer_y*highLNI",
                      reg.data = OutMigrationData)

reg_out <- str_replace(reg_out, "SHI25:farmer\\\\_y:highLNI", 
                       "\\\\multirow[t]{3}{*}{\\\\makecell[l]{\\\\\\\\ \\\\\\\\ \\\\\\\\Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer $\\\\times$ \\\\\\\\ \\\\qquad High Communal \\\\\\\\ \\\\qquad Identification}}")
reg_out <- str_replace(reg_out, "SHI25:farmer\\\\_y", 
                       "\\\\multirow[t]{2}{*}{\\\\makecell[l]{\\\\\\\\ Soil Heterogeneity $\\\\times$ \\\\\\\\ \\\\quad Farmer}}")
reg_out <- str_replace(reg_out, "SHI25", "Soil Heterogeneity")
reg_out <- str_replace(reg_out, "farmer\\\\_y", "Farmer")

outfile <- paste0(AppendixTablesPath, "/table_D_21.tex")
write(c(reg_out[c(3:11,18:20)], "  & & & & & & \\\\ ", 
        "  & & & & & & \\\\ ", "\\hline \\\\[-1.8ex]  \\\\[-1.8ex] ", 
        reg_out[21:22]), outfile)


#--- Foot ---
footLines <- 
  paste(
    c(
      paste("State $\\times$ Year Fixed Effects", paste(rep(" ", 3), collapse = " & "), paste(rep("\\checkmark", 5), collapse = " & "), "\\\\"),
      paste("Geoclimatic Controls", paste(rep(" ", 4), collapse = " & "), paste(rep("\\checkmark", 4), collapse = " & "), "\\\\"),
      paste("Smooth Location Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Agricultural Suitability Controls", paste(rep(" ", 5), collapse = " & "), paste(rep("\\checkmark", 3), collapse = " & "), "\\\\"),
      paste("Higher Order Controls", paste(rep(" ", 6), collapse = " & "), paste(rep("\\checkmark", 2), collapse = " & "), "\\\\"), 
      paste("Individual Controls", paste(rep(" ", 7), collapse = " & "), paste(rep("\\checkmark", 1), collapse = " & "), "\\\\")), 
    collapse = "\n") 

outfile <- paste0(AppendixTablesPath, "/table_D_21_foot.tex")
write(unlist(footLines), outfile)

# Clear -------------------------------------------------------------------

rm(list = ls())
gc()